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The  Bifurcation  Interpreter: 

A  step  towards  the  automatic  analysb 
of  dynamical  systems 

Haxold  Abelson 


Abstract 

The  Bifurcation  Interpreter  is  a  computer  program  that  autonomously 
explores  the  steady-state  orbits  of  one-parameter  families  of  periodi<'ally- 
driven  oscillators.  To  report  its  findings,  the  Interpreter  generates 
schematic  diagrams  and  English  text  descriptions  similar  to  those  ap¬ 
pearing  in  the  science  and  engineering  research  literature.  Given  a 
system  of  equations  input,  the  Interpreter  uses  symbolic  algebra 
to  automatically  generate  numerical  procedures  that  simulate  the  sys¬ 
tem.  The  Interpreter  incorporates  knowledge  about  dynamical  sys¬ 
tems  theory,  which  it  uses  to  guide  the  simulations,  to  interpret  the 
results,  and  to  minimize  the  effects  of  numerical  error. 

t. 

Most  of  the  effort  in  scientific  computing  concerns  programs  that  pro¬ 
duce  numerical  output,  although  there  is  also  much  interest  in  developing 
graphical  rendering  techniques  that  help  people  to  visualize  this  output  [11]. 
But  whether  a  numerical  experiment  generates  numbers  or  pictures,  there  re¬ 
mains  the  task  of  interpreting  the  results — to  distill  the  numerical  output  into 
high-level,  often  quaditative  descriptions  that  can  be  summarized,  reasoned 
about,  and  used  to  guide  new  experiments.  This  paper  illustrates  how  such 
interpretations  can  be  created  automat icadly,  with  appropriate  combinations 
of  numerical  and  symbolic  processing. 

The  Bifurcation  Interpreter  is  a  computer  program  that  autonomously 
explores  the  steady-state  orbits  of  one- parameter  families  of  periodically- 
driven  oscillators.  To  report  its  findings,  the  Interpreter  generates  schematic 
diagrams  and  English  text  descriptions  similar  to  those  appearing  in  the 
science  and  engineering  research  literature. 

Section  one  of  this  paper  illustrates  reports  generated  by  the  Interpreter 
and  compares  these  with  similar  reports  published  by  dynamicists  to  de¬ 
scribe  their  own  investigations  of  nonhnear  systems.  Section  two  sununa- 
rizes  the  range  of  phenomena  that  the  Bifurcation  Interpreter  can  obsCTve, 
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and  describes  the  internal  data  structures  it  maintains  to  keep  track  of  its 
observations.  Section  three  is  a  detadled  description  of  the  methods  used  by 
the  program.  The  paper  ends  by  discussing  the  limitations  of  the  present 
implementation  and  describes  directions  for  improvement. 

1  Qualitative  behavior  of  dynamical  systems 

A  periodically-driven  dissipative  nonlinear  oscillator  typically  admits  a  finite 
number  of  periodic  orbits,  where  the  period  of  each  orbit  is  an  integer  multiple 
of  the  drive  ^  TV'S!;  multiple  is  called  the  ordci  of  ifUC-  Oiblt.  IaA  a 

parametrized  family  of  oscillators,  smooth  variations  in  the  parameters  give 
rise  to  families  of  orbits  that  also  vary  smoothly,  except  at  certain  critical 
values  of  the  parameters  at  which  there  occur  bifureationSy  or  discontinuous 
changes  in  orbital  behavior.  Deprading  on  the  t)rpe  of  the  bifurcation,  a 
family  of  orbits  may  vanish,  change  order,  or  merge  with  other  families,  or 
new  families  of  orbits  may  be  bom.  Families  that  meet  at  a  bifurcation  are 
said  to  be  in  the  same  class. 

The  geometric  theory  of  dynamical  systems  focuses  on  the  evolution  of 
steady-state  orbits  through  bifurcations  as  a  framework  for  capturing  the 
qualitative  behavior  of  dynamical  systems.  A  major  achievement  of  the  the¬ 
ory  has  been  to  show  that  for  one-parameter  families,  at  least,  the  bifurca¬ 
tions  typically  encountered  can  be  classified  into  a  small  number  of  recog¬ 
nized  types;  the  type  of  the  bifurcation  determines,  up  to  homeomorphism, 
the  behavior  of  the  faunily  near  the  bifurcation  point.* 

1.1  A  sample  report  by  the  Interpreter 

Given  a  one-parameter  family  of  Deriodically-driven  oscillators,  the  Bifurca¬ 
tion  Interpreter  identifies  fan.  .;>es  ..f  stable  periodic  orbits  and  classifies  the 

Mn  addition,  nonlinear  oscillators  tuay  have  non-periodic  steady-state  orbits.  These 
include  luasiperiodic  orbits,  which  have  discrete-frequency  spectra,  but  not  at  rational 
multiples  of  the  drive  frequency;  and  ckaotic  orbits,  which,  loosely  speaking,  are  steady- 
state  orbits  that  are  neither  periodic  nor  quasi-periodic. 

’Holmes  and  Guckenheimer  [6]  present  the  mathematical  foundations  of  bifurcation 
theory.  Thompson  and  Stewart  [19]  provide  an  introduction  to  dynamical  systems  theory 
for  scientists  and  engineers. 
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bifurcations  through  which  they  evolve.  The  Interpreter’s  report  catalogues 
the  classes  of  orbits  and  the  types  of  the  bifurcations. 

As  iin  example,  the  following  form  of  Buffing’s  equation 

X  +  kx  +  =  p  cos  t 

describes  a  one-parameter  family  of  nonlinear  oscillators  for  each  fixed  value 
of  k.  This  family  for  A:  =  0.1  was  presented  to  the  Bifurcation  Interpreter  for 
analysis  by  specifying  the  equation  (as  a  system  of  two  first-order  equations), 
a  varying  parameter  and  parametci  interval,  and  a  domain  in  state-space 
outside  of  which  orbits  should  be  considered  to  be  unbounded: 

•  system: 

Xi  =  X2t  ^2  =  pcos  t  —  kx2  —  (1) 

•  fixed  parameters:  ^  =  0.1 

•  varying  parameter:  p,  1  to  25 

•  bounds  on  state  variables:  ij,  —5  to  5;  xj,  —10  to  10 

The  Bifurcation  Interpreter’s  output  is  an  automatically-generated  tex¬ 
tual  report  that  presents  the  results  of  the  analysis  using  the  technical  terms 
of  the  dynamics  literature.  (The  meemings  of  these  terms  are  explained  in 
section  2  below.)  Here  is  an  excerpt: 

The  system  was  explored  for  values  of  p  between  1  and  25,  and  10 
classes  of  stable  periodic  orbits  were  identified. 

Class  A  is  already  present  at  the  start  of  the  parameter  range  p  =  1 
with  a  family  of  order-i  orbits  Ao.  Near  p  =  2.287,  there  is  a  super¬ 
critical  pitchfork  bifurcation,  and  Ao  splits  into  symmetric  families 
Aifl  and  Ai,i,  each  of  order  1.  Ai,o  vanishes  at  a  fold  bifurcation  near 
p  =  3.567.  Ai,i  vanishes  similarly. 

Class  B  appears  around  p  =  3.085  with  a  family  of  order- 1  orbits 
Bo  arising  from  a  fold  bifurcation.  As  the  parameter  p  increases,  Bq 
undergoes  a  period-doubling  cascade,  reaching  order  2  near  p  =  4.876, 
and  order  4  near  p  =  5.441.  Although  the  cascade  was  not  traced  past 
the  order  4  orbit,  there  is  apparently  another  period-doubling  near 
p  =  5.52,  and  a  chaotic  orbit  was  observed  at  p  =  5.688. 
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Class  D  appears  around  p  =  5.642  with  a  family  of  order-3  orbits 
Do  arising  from  a  fold  bifurcation.  Near  p  *  6.691,  there  is  a  supercrit¬ 
ical  pitchfork  bifurcation,  and  Dq  splits  into  symmetric  families  Di^ 
and  each  of  order  3.  Near  p  =  7.992,  there  is  a  sup^critical- 
pitchfork  bifurcation  where  Dxjo  and  Di^i  merge  to  form  a  family  Dj 
of  order  3.  Near  p  =  9.677,  there  is  a  transcritical  bifurcation  where 
D]  vanishes  by  exchangpng  stability  with  family  D3  of  order  3.  D3 
vanishes  at  a  bifurcation  of  undetermined  type  near  p  =  9.858. 

Class  J  appears  around  p  =  23.96  as  a  family  of  order-5  orbits 
Jo  arising  from  a  fold  bifurcation.  Jo  is  present  at  the  end  of  the 
parameter  range  at  p  =  25. 

In  addition  to  the  text  report,  tne  Interpreter  can  generate  diagrams,  such 
as  the  one  in  figure  1,  that  present  its  findings  in  schematic  form.^  As  shown 
in  the  figure,  the  horizontal  axis  represents  the  parameter  range,  and  tic- 
marks  indicate  the  values  where  bifurcations  occur.  The  classes  are  labelled 
and  arrayed  verticaUy  to  show  the  parameter  extent  over  which  orbits  in 
each  class  occur.  For  example,  for  values  of  p  between  pa  «  3.1  and  «  3.6, 
the  system  has  four  coexisting  steady- state  orbits  (two  in  class  A,  one  in  B, 
and  one  in  C);  each  system  trajectory  will  evolve  to  one  of  these  four  orbits, 
depending  upon  the  initial  conditions.  Each  bifurcation  type  is  indicated 
by  a  standaird  symbol — the  bracket  indicates  a  fold,  the  “Y”  a  pitchfork, 
and  the  question  mark  within  a  circle  a  bifurcation  that  the  Interpreter  was 
unable  to  identify.  The  blank  rectangles  denote  parameter  intervals  over 
which  the  diagram  should  be  expanded  in  order  to  reveal  further  details. 
Figure  2  shows  this  expansion  for  the  indicated  regions  in  classes  B  and  C, 
and  figure  3  shows  details  of  cl2is8  D. 

Observe  that  the  Interpreter’s  textual  report  includes  information  that  is 
not  reflected  in  the  diagrams.  For  instance,  the  Interpreter  has  recognized 
the  succession  of  bifurcations  in  class  B  as  a,  period-doubling  cascade,  and 
noted  that  the  existence  of  the  cascade  is  consistent  with  its  observation  of 
a  chaotic  orbit  at  a  higher  value  of  the  parameter. 

^The  graphical-output  routines  used  in  the  Interpreter  program  were  implemented  by 
Ognen  Nastov,  under  the  auspices  of  the  MIT  Undergraduate  Research  Opportunities 
Program. 
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f(Jd  supercritical  pitchfork  unidentified  bifurcation 
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supercritical  flip  expand  for  further  detail 


Figure  1:  This  diagram,  generated  automatically  by  the  Bifurcation  Interpreter,  shows 
the  analysis  of  the  system  described  by  equation  (1)  for  k  =  0.1.  The  Interpreter  has 
traced  the  evolution  of  10  classes  of  families  of  periodic  orbits  and  their  bifurcations.  Each 
bifurcation  type  is  identified  by  a  standard  symbol — the  key  indicates  the  Interpreter’s 
symbols  for  fold  bifurcations,  supercritical-pitchfork  bifurcations,  and  for  bifurcations  that 
the  Interpreter  has  failed  to  identify.  The  blank  rectangles  indicate  parameter  intervals 
over  which  the  figure  should  be  expanded  to  reveal  mote  detail,  as  shown  in  figures  2 
and  3.  The  p  values  along  the  horizontal  axis  indicate  the  parameter  values  at  which  the 
bifurcations  occur. 
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Fig^ure  2:  Expanding  the  interval  from  p  =  3  to  p  =  5.7  showi  detaila  of  claaaea  B  and 
C  from  figure  1.  Each  claa.  begins  as  order  1,  then  doubles  to  order  2  vU  aupercritical- 
flip  bifurcations.  Expanding  the  interval  indicated  by  the  rectangles  would  show  further 

doubUng.  The  Interpreter  recognizes  this  progreaaion  as  a  period-doubling  cascade  and 
repiorts  it  as  such.  ’ 


Figure  3:  Expanding  the  interval  from  p  =  5.6  to  p  =  10  shows  details  of  the  class  D 
of  order-3  orLits.  Starting  as  a  single  family,  the  clawi  sphts  into  a  pair  near  p  =  6.7  and 
these  merge  again  near  p  =  8. 
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1.2  Comparison  with  manually-produced  reports 

The  diagrams  and  the  text  generated  by  the  Bifurcation  Interpreter  are  sim¬ 
ilar  to  accounts  that  scientists  and  engineers  have  published  to  document 
their  investigations  of  nonlinear  systems.  One  example,  shown  in  figure  4,  is 
a  diagram  presented  by  V.  Franceschini  [5]  to  report  the  results  of  a  series 
of  simulation  studies  of  fluid  flow,  modeled  using  a  one-parameter  family  of 
nonlinear  equations  parametrized  by  the  Reynolds  number.  Franceschini  de¬ 
scribes  his  diagram  as  a  “graphical  summairy  of  the  phenomenology”  of  the 
model,  and  compares  diagrams  for  different  models  to  determine  whether 
they  exhibit  the  same  behavior."*  Although  the  detailed  layout  of  the  dia¬ 
gram  is  different  from  those  generated  by  the  Bifurcation  Interpreter,  the 
information  presented  is  much  the  same.®  Incidentally,  since  the  Interpreter 
generates  its  diagrams  from  a  symbolic  data  structure  that  it  maintains  in¬ 
ternally  (described  in  section  2.2  below),  it  should  not  be  difficult  to  extend 
the  Interpreter  to  automatically  contrast  the  qualitative  behaviors  of  two 
different  systems  as  does  Franceschini  in  his  paper. 

The  Interpreter’s  reports  are  based  upon  numerical  experiments  con¬ 
ducted  at  finite  resolution,  so  it  is  possible  that  the  program  may  fail  to 
observe  orbits  that  exist  only  over  small  parameter  intervals  or  small  regions 
of  state  space.  On  the  other  hand,  experiments  conducted  by  people  suffer 
the  same  criticism.  Indeed,  one  can  check  the  qu2dity  of  the  Interpreter’s 
observations  of  the  system  in  equation  (1)  by  comparing  them  with  an  often- 
cited  study  published  by  Y.  Ueda  [20],  who  has  catalogued  the  steady-state 
solutions  of  this  system  over  the  range  0  <  k  <  0.8,  0  <  p  <  26  by  means  of 

*The  Navier-Stokes  equations  for  fluid  flow  are  a  system  of  partial  differential  equa- 
Expaiiding  liiC  solution  a  Fourier  series  and  truncating  the  expansion  at  a  fixed 
order  reduces  this  to  a  system  of  ordinary  differential  equations.  One  approach  to  study¬ 
ing  the  Navier-Stokes  equations  is  to  study  these  finite-mode  approximations  in  the  hope 
that  their  solutions  will  have  the  same  qualitative  properties  as  those  of  the  full  equations. 
Franceschini ’s  experiments  [5]  test  the  validity  of  this  approach  by  comparing  the  quali¬ 
tative  properties  of  the  8-mode  and  9-mode  models.  He  finds  that  the  two  models  have 
almost  completely  different  phenouiecologies,  althcu^h  ran  p-biKit  turbulence  at  low 
Reynolds  numbers. 

®The  most  significant  difference  between  the  two  types  of  diagrams  is  that  figure  4 
indicates  unstable  periodic  orbits  (the  dotted  line  marked  beginning  at  the  lower  left  of 
the  figure)  and  shows  regions  of  chaotic  behavior  (at  the  right-hand  sides  of  the  main  figure 
and  of  insert  E).  The  Bifurcation  Interpreter  does  not  currently  record  such  information, 
but  could  be  extended  to  do  so. 


Figure  4:  This  diagram,  reproduced  from  [5],  ia  a  manually-produced  report  aimilar  to 
the  diagrams  that  the  Bifurcation  Interpreter  generates  automatically. 


extensive  simulations  with  analog  and  digital  computers.  In  a  series  of  tests 
with  various  values  of  k,  the  Interpreter  identified  all  the  orbits  in  Ueda’s  cat¬ 
alogue  and  sometimes  observed  orbits  that  are  missing  from  the  catalogue. 
For  instance,  the  )rder-5  family  J  shown  above  in  the  Interpreter’s  report 
for  it  =  0.1  is  absent  from  Ueda’s  catalogue.  Also,  although  Ueda  describes 
the  order-3  orbit  in  cIms  D,  he  does  not  note  the  split  into  distinct  families 
between  p  =  6.7  and  p  =  8  that  is  indicated  in  figure  3.  Figures  5  and  6 
illustrate  these  missing  orbits. 


2  What  the  Interpreter  can  recognize 

The  Interpreter  analyzes  one-parameter  families  of  dynamical  systems.  These 
can  be  discrete  dynamical  systems,  formed  by  iterating  diffeomornhisms  fp 
(p  is  the  family  parameter);  or  they  can  be  systems  of  ordinary  differential 
equations  ir  =  Hp{x,  t),  where  each  ffp  is  periodic  in  t  with  period  The 


®  At  present,  the  Interpreter  has  been  fully  implemented  only  for  systems  with  a  two- 
dimensional  state  space  x  =  (xi,  X]),  although  most  of  the  underlying  algorithms  will  work 


Figure  5:  Shown  here  and  in  figure  6  are  classes  of  steady-state  orbits  automatically 
reported  by  the  Bifurcation  Interpreter  for  the  system  in  equation  (1)  that  are  missing 
from  a  well-known  published  catalogue  of  orbits  for  this  system  [20].  Plotted  here  are  x{t) 
and  the  drive  for  an  order-5  orbit  from  class  J  of  figure  1  (fc  =  0.1,p  =  24.25). 

periodic  continuous  case  reduces  to  the  discrete  caise  by  the  standard  trick 
of  considering  fp  to  be  the  period  map,  which  maps  each  state  x  to  the  point 
t  =  tff  on  the  trajectory  of  Hp  that  starts  from  x  at  t  =  0.  Since  Hp 
is  periodic,  the  successive  iterates  are  the  points  on  the  trajectory 

at  times  that  axe  multiples  of  tfj-  Thus,  a  trajectory  that  is  periodic  with 
period  ntn  corresponds  to  periodic  point  of  order-n  point  of  the  period  map, 
that  is,  a  point  i  such  that  /^"^x)  =  x. 

Given  a  system  to  analyze,  a  parameter  interval  pi  <  p  <  Ph,  and  a 
bounded  domain  in  state  space,  the  Interpreter  attempts  to  determine,  for 
each  p  in  the  parameter  interval,  the  periodic  points  of  fp  that  lie  within  the 
specified  domain.  It  does  this  by  locating  periodic  points  for  a  few  v2Llues  of 
p,  and  tracing  the  family  of  periodic  points  that  evolves  as  p  varies.  A  family 
may  persist  until  the  end  of  the  parameter  interval;  or  a  family  may  vanish 
because  the  periodic  point  goes  out  of  bounds  (leaves  the  designated  state- 
space  domain);  or  the  family  may  vanish  when  the  periodic  point  undergoes 
a  bifurcation. 

in  any  number  of  dimensions.  See  the  discussion  in  section  4. 
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Figure  6:  This  is  the  phase  portrait  for  one  of  two  distiDCt  order-S  orbits  in  class  D 
(figure  1)  that  coexist  at  the  same  parameter  values  (ib  =  0.1, p  =  7.45),  a  phenomenon 
reported  by  the  Bifurcation  Interpreter,  but  not  noted  in  the  catalogue  published  in  [20]. 
The  phase  portrait  of  the  second  orbit  is  the  reflection  of  this  one  under  the  symmetry 
{x,  x)  •-+  {—X,  —x).  The  crosses  drawn  on  the  orbit  indicate  the  Poincari  poinU,  at  which 
<  is  a  multiple  of  the  period. 


2.1  Types  of  bifurcations 

A  bifurcabion  of  a  periodic  point  x  of  order  n  can  be  recognized  computa- 
tionadly  by  the  fact  that,  as  p  vzo-ies,  an  eigenvalue  of  D the  Jaw:obian 
derivative  of  fp  at  i,  crosses  the  unit  circle  in  the  complex  planed  When 
the  Interpreter  encounter  \  bifurcation,  it  attempts  to  classify  the  bifurca¬ 
tion,  based  upon  the  eigenvalues  of  jD/^”^(x)  and  the  pattern  of  stable  and 
unstable  periodic  points  meeting  at  x. 

Here  are  the  types  of  bifurcations  that  the  Interpreter  recognizes;* 

•  fold  bifurcation:  For  values  of  p  to  one  side  of  the  bifurcation  value  a 
stable  and  an  unstable  periodic  point  (of  the  same  order)  coexist.  These 
two  periodic  points  meet  as  p  approaches  the  bifurcation  value;  on  the 
other  side  of  the  bifurcation,  both  periodic  points  have  vanished.®  An 
eigenvalue  of  Df^^\x)  crosses  the  unit  circle  at  1. 

•  flip  bifurcations: 

—  supercritical  flip:  A  stable  periodic  point  of  order  n  sphts  to  form 
an  unstable  periodic  point  of  order  n  and  a  stable  periodic  point 
of  order  2n. 

-  subcritical  flip:  This  is  equivalent  to  a  supercritical  flip  of  /”*,  the 
inverse  of  /;  an  unstable  periodic  point  of  order  2n  merges  with  a 
stable  periodic  point  of  order  n  to  form  an  unstable  periodic  point 
of  order  n. 

An  eigenvaJue  of  Z?/^"^(x)  crosses  the  unit  circle  at  —1. 

•  Niemark  bifurcations: 

^So  long  as  no  eigenvalue  has  modulus  1,  the  periodic  point  x  is  hyperbolic  and  thus 
will  not  change  type  with  small  changes  in  p.  See  [6]  for  details. 

*In  referring  to  the  bifurcation  types,  we  have  adopted  the  terminology  of  [19].  Various 
authors  use  different,  and  even  incompatible  terminology.  For  example,  the  flip  is  also 
called  a  cusp  or  a  pitchfork. 

^Depending  on  the  direction  from  which  the  fold  is  encountered,  it  wiU  be  observed 
either  as  a  stable  periodic  point  and  an  unstable  periodic  point  meeting  and  annihilating 
each  other,  or  alternatively  as  the  birth  of  a  stable-unstable  pair  of  periodic  points.  Since 
the  Interpreter  ordinarily  encounters  a  bifurcation  by  tracking  a  stable  point  into  it,  we 
describe  the  bifurcations  from  that  orientation. 
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--  supercritical  Niemark:  A  stable  periodic  point  vanishes,  forming 
an  unstable  periodic  point  and  a  stable  limit  cycle. 

-  subcritical  Niemark:  A  stable  periodic  point  meets  an  unstable 
limit  cycle,  resulting  in  an  unstable  periodic  point. 

An  eigenvalue  of  crosses  the  unit  circle  at  a  value  other  than 

1  or  —1. 

•  transcritical  bifurcation:  A  stable  periodic  point  and  an  unstable  pe¬ 
riodic  point  meet  and  “exchange  stabilities":  on  the  other  side  of  the 
bifurcation,  the  (extrapolated)  stable  point  is  now  unstable  and  the 
unstable  point  is  stable.  An  eigenvalue  of  Df!^\x)  crosses  the  unit 
circle  at  1. 

•  pitchfork  bifurcations: 

-  supercritical  pitchfork:  A  stable  periodic  point  splits  to  form  two 
stable  periodic  points  and  an  unstable  periodic  point,  all  of  the 
same  order. 

-  subcritical  pitchfork:  An  stable  periodic  point  merges  with  two 
unstable  periodic  points,  resulting  in  an  unstable  periodic  point. 

An  eigenvalue  of  Df^\x)  crosses  the  unit  circle  at  1. 

These  bifurcations  are  summarized  in  figure  7  by  means  of  bifurcation 
diagrams,  of  the  kind  commonly  employed  in  the  dynamics  literature.  These 
diagrams  show,  near  each  bifurcation,  the  paths  of  the  periodic  points  as 
p  varies,  projected  onto  a  one-dimensionaJ  subspace  of  phase  space  (two- 
dimensional  subspace  for  the  Niemark  bifurcations).  The  parabolic  paths 
drawn  for  most  of  the  bifurcations  are  more  than  simply  schematic.  It  can 
be  shown  (see  [19])  that,  for  values  of  p  very  near  the  bifurcation  value,  the 
periodic  points  approach  the  bifurcations  along  such  pairabolic  paths.  The 
methods  by  which  the  Interpreter  recognizes  bifurcations  maike  use  of  this 
fact  (section  3.4). 

The  bifurcations  listed  here  are  the  typirjil  ones  encountered  in  one- 
parameter  families  of  maps,  in  the  following  sense:  It  can  be  shown  that 
the  only  bifurcations  that  appear  in  generic  one-parameter  families  are  the 
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Supercritical  Subcritical  Supercritical 

Pitchfork  Pitchfork  Niemark 


Subcritical 
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Figure  7:  This  ia  the  set  of  bifurcations  types  recognized  by  the  Bifurcation  Interpreter. 
The  figure  shows  commonly  used  “bifurcation  diagrams”  that  describe  the  pattern  of  stable 
and  unstable  periodic  points  near  each  bifurcation,  in  which  the  horizontal  axis  represents 
the  parameter  direction  and  the  vertical  axis  (verical  plane  in  the  case  of  the  Niemark 
bifurcations)  represents  the  state-space.  The  curve  indicate  the  loci  of  the  period!  points 
as  the  parameter  varies.  Open  circles  indicate  the  loci  of  unstable  points  and  solid  circles 
indicate  stable  points. 
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fold,  the  flip,  and  the  Niemark  bifurcations.  In  addition,  the  transcriticai  bi¬ 
furcation  is  generic  for  one- parameter  families  of  maps  that  are  constrained 
to  have  a  common  fixed  point  for  all  fp;  and  the  pitchfork  bifurcations  are 
generic  in  families  of  maps  that  are  invariant  under  a  symmetry.  See  [6] 
and  [19]  for  more  details. 

In  addition  to  recognizing  individual  bifurcations,  the  Interpreter  also  no¬ 
tices  typical  patterns  of  bifurcations.  One  such  pattern  is  the  period- doubling 
cascade  [4],  in  which  an  infinite  sequence  of  supercritical-flip  bifurcations  oc¬ 
curring  at  a  converging  sequence  of  values  of  p  produces  periodic  points  of 
orders  n,  2n,  4n, . . . ,  followed  by  chaotic  orbits  at  values  of  p  beyond  the  ac¬ 
cumulation  point.  A  second  pattern  is  the  symmetric  pair,  produced  when 
an  orbit  that  is  invziriant  under  a  symmetry  of  the  system  loses  its  synametry 
and  splits,  at  a  supercritical-pitchfork  bifurcation,  into  two  orbits,  each  of 
which  is  transformed  into  the  other  by  the  symmetry.^®  Continued  evolution 
of  these  orbits  produces  symmetric  pairs  of  families  and  bifurcations. 

2,2  The  family  report 

Although  the  text  and  diagrams  illustrated  in  section  1.1  are  the  must  con¬ 
spicuous  part  of  the  Bifurcation  Interpreter’s  output,  the  bulk  of  the  program 
is  concerned  with  conducting  numerical  experiments  and  distilling  their  re¬ 
sults  into  a  data  structure  called  a  family  report.  The  family  report  is  the 
basis  for  generating  the  text  and  diagrams.  It  can  also  be  examined  by  other 
programs  to  answer  questions  about  dynamical  systems,  for  example  to  in¬ 
dicate  the  ranges  over  which  multiple  families  exist,  or  to  list  the  types  of 
bifurcations  encountered. 

The  family  report  lists  the  classes  of  families  in  order  of  appearance  with 
increasing  p.  For  instance,  the  analysis  of  DufiBng  system  shown  in  figure  1 
has  10  classes.  The  first  of  these,  class  A,  contains  three  families,  two  of 
which  form  a  symmetric  pair: 

(class :1 
(■sabsrs 

(  #Cfa«lly:A.O]  (symstric-palr  #Cla»lly:A.1.03  ittaallyiA.l.!])))) 

^°For  example,  the  Duffing  system  described  in  section  1.1  has  the  form  s  = 
where  H(—x,t  -t-  r)  =  -ff(x,i).  Once  can  show  that  this  entails  a  symmetry:  for  every 
closed  orbit  5,  the  set  -5  =  {-*  |  at  6  5}  is  also  a  closed  orbit.  Thus,  either  5  itself  is 
invariant  under  the  symmetry,  or  is  one  of  a  symmetric  pair  of  closed  orbits,  as  in  figure  6. 
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Each  family  is  represented  by  a  data  structure  that  contains  all  the  infor¬ 
mation  collected  about  the  family,  such  eis  the  name  and  order  of  the  family, 
the  bifurcations  at  which  the  family  appears  and  vanishes,  and  the  instances 
of  the  periodic  point  that  was  tracked  to  form  the  family.  Here  is  the  entry 
for  Ai,o,  the  second  family  in  class  A.  This  entry  includes  the  information 
that  Ai,o  is  part  of  a  symmetric  pair,  paired  with  faunily 
(family:!. 1.0 
(ordar  1) 

(appears  #  [bifurcation :bif.l]) 

(vanishes  #[bifurcation:bif .6]) 

(paired-vith  (#[fajaily:A.l.l])) 

(instances  69  #[list  of  instances  ...])) 

The  list  of  instances  describes  the  individual  periodic  points.  Here  is  the 
first  point  in  the  list  for  Ai,o 
(instance 

(parasieter-value  2.2867) 

(type  periodic) 

(orbit  «(1.8861  .36982)) 

(jacobiaa  #(#(-.29979  -9.667ee-2)  #(11.327  1.8317))) 

(eigenvalues  .99666  .63628) 

(sensitivity  #(-7.9680  111.496))) 

This  is  a  periodic  point  at  p  =  2.2867  and  i  =  (1.8851, 0.3998).  Also  recorded 
are  the  Jacobian  amd  the  eigenvalues  of  the  period  map  at  this  point,  and 
the  sensitivity  derivative  (i.e.  the  derivative  with  respect  to  p  (vector- valued) 
function  p  f{p,x)),  which  is  used  in  tracking  periodic  points  as  explained 
in  section  3.3. 

In  addition  to  the  list  of  classes,  the  Interpreter  generates  a  list  of  the 
bifurcations  encountered.  Here  is  the  bifurcation  at  which  appears: 

(bif urcat ion ; bif . 1 
(typa  snpsrcrltical-pitchfork) 

(paramstsr-rangs  (2.2864  2.2867)) 

(famillss-vanishing  (# [family: A. 0] )) 

(familiss-appaaring  (# [family: A. 1.0]  # [family: A. 1 . 1] )) 

(symmatrlc-pair  (split  to  form  symmatric  pair)) 

(dlractlon  paramatar-incraasing) 

(coaoiaats  ((infarrad  from  marging  familias)))) 

This  is  a  supercritical-pitchfork  bifurcation  located  between  p  =  2.2864 
and  p  =  2.2867,  where  Ao  vanishes  2ind  Ai^  and  Ai,i  appear.  The  vamishing 
family  splits  to  form  a  symmetric  pair  of  families  with  a  common  head  (in- 
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stance  at  lowest  value  of  p).  The  bifurcation  is  oriented  in  the  “parameter 
increasing”  direction  (that  is,  the  split  happens  as  p  increases).  Finally,  a 
comment  records  that,  in  this  case,  the  Interpreter  discovered  the  bifurcation 
by  noticing  that  two  families  merged  (see  figure  10b). 

Finally,  the  Interpreter  maintains  a  list  of  the  parameter  values  at  which 
it  found  orbits  that  it  was  not  able  to  classify.  Depending  on  the  nearby 
bifurcations,  these  might  be  interpreted  as  chaotic  orbits.  Given  the  infor¬ 
mation  in  the  family  report,  it  should  be  clear  how  both  the  text  and  the 
diagrams  in  section  1.1  can  be  generated. 

3  How  the  Interpreter  works 

Figure  8  shows  the  major  modules  in  the  Bifurcation  Interpreter  and  the 
information  that  flows  from  one  module  to  the  next.  B^inning  with  equar 
tions  that  describe  a  dynamical  system,  and  a  designated  parameter  range 
to  explore,  the  Interpreter  generates  numerical  procedures  that  compute  the 
transformation  from  state  to  state,  the  Jacobian  of  the  transformation,  and 
the  sensitivity  (the  derivative  of  the  transformation  with  respect  to  the  vary¬ 
ing  parameter).  With  the  aid  of  these  procedures,  the  Interpreter  locates 
stable  periodic  points  at  various  parameter  values.  Each  periodic  point  is 
tracked  as  the  parameter  varies,  producing  a  family  that  generaUy  ends  in  a 
bifurcation.  The  Interpreter  attempts  to  classify  the  bifurcation  based  upon 
information  obtained  in  the  tracking  process.  Typically,  classifying  the  bi¬ 
furcation  will  produce  new  stable  periodic  points  that  must  themselves  be 
tracked. 

When  all  known  periodic  points  have  been  extended  to  complete  fami¬ 
lies,  the  Interpreter  reexamines  the  bifurcations  in  the  light  of  more  global 
information,  and  corrects  its  previous  classifications  if  necessary.  Next,  bifur¬ 
cations  and  families  are  examined  for  symmetric  pairs  and  period-doubling 
cascades.  The  Interpreter  consolidates  families  into  classes,  assigns  names  to 
the  families  and  bifurcations,  and  organizes  the  data  into  a  family  report. 
The  family  report  is  the  basis  for  generating  text  and  diagrams  as  exhibited 
in  section  1.1. 

We  now  consider  each  of  these  steps  in  detail. 


Figxire  8:  The  block  diagram  shows  the  major  components  of  the  Bifurcation  Interpreter, 
begiiming  with  a  system  specified  by  algebraic  equations,  and  ending  with  a  family  report, 
from  which  text  and  diagrams  are  generated.  The  arrows  between  blocks  indicate  the 
principal  output  of  each  module,  which  is  the  input  to  the  next  module. 


18 


3.1  Compiling  procedures  for  Newton’s  method 

Most  of  the  Interpreter’s  algorithms  rely  upon  the  ability  to  locate  fixed 
point®  This  is  accomplished  with  the  £ud  of  Newton’s  method.  Applying 
Newton’s  method  to  find  fixed  points  of  a  map  /  requires  evaluating  both  / 
and  the  Jacobian  derivative  Df  at  various  points  x  in  the  state  space.  For  a 
discrete  dynamical  system,  where  /  is  specified  explicitly  by  algebraic  equa¬ 
tions,  the  Interpreter  differentiates  the  equations  symbolically  to  produce 
expressions  for  the  components  of  the  Jacobian  matrix  [Df{x)]ij  =  dpjdxj. 
These  expressions  are  compiled  into  a  procedure  that  evaluates  the  Jacobian, 
given  a  numerical  value  for  x. 

For  a  system  of  ordinary  differential  equations  x  =  H{x,t)  with  H{x,t) 
periodic  in  t,  /  is  not  given  by  an  explicit  algebraic  equation.  Rather,  / 
is  the  period  map,  evaluated  by  integrating  H  along  the  trajectory  start¬ 
ing  from  X.  The  components  of  the  (transpose  of  the)  Jacobian  matrix, 
<.7  =  [0/(^)15  =  [DS{x)\i^  are  computed  by  integrating,  along  the  ssune 
trajectory,  the  associated  variational  system 

dt  ~^dxu 

with  initial  conditions  ^,j(0)  =  0  for  i  ^  j  and  ^ti(O)  =  1  for  »  =  j.  The 
Interpreter  symbolically  differentiates  the  expressions  for  H{x,t)  to  obtain 
expressions  for  the  dH'  jdxj  and  compiles  these  into  a  procedure  for  com¬ 
puting  the  Jacobian  by  numeric2J  integration. 

For  example,  figure  9  shows  the  DufiBng  system  of  equation  (1)  as  it  is 
actually  input  to  the  Interpreter,  with  a  specification  of  the  state  variables, 
the  parameter,  and  equations  for  the  derivatives  of  the  state  variables.^^ 
The  Interpreter  produces  from  this  a  system-derivative  generator.  This  is  a 
procedure  which,  given  values  for  the  parameters  k  and  p,  returns  a  system 
derivative:  a  procedure  which,  given  a  state  vector  (t,  Xj,  Z2),  returns  a  vector 
whose  components  are  the  derivatives  with  respect  to  i  of  the  state-vector 
components.  Figure  9  also  shows  the  procedure  generated  for  evolving  the 
associated  variational  system.  Given  values  for  k  and  p,  this  produces  a 

^^The  input  language  used  here,  phrased  in  terms  of  constraints,  is  more  general  than 
is  required  for  the  purposes  of  the  Interpreter.  The  same  language,  and  the  techniques 
for  generating  system  derivatives  and  variational  equations,  were  used  by  Abelson  and 
Sussman  to  support  an  electrical-circuit  an^ysis  program  as  described  in  [1]. 
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vajiational-system  derivative,  which  operates  on  a  variational  state  vector 
with  components  f,  ii,  and  the  returning  the  vector  whose  components 
are  the  time  derivatives  of  t,  xi,  12,  together  with  the  components  of  the 
Jacobian,  computed  as 

— —  k6\2 
S22 

— 3Xi^31  —  ^<^23 

The  Interpreter  generated  the  expressions  for  these  derivatives  by  symbolic 
differentiation. 

The  system-derivative  procedures  generated  by  the  Interpreter  c^ul  be 
combined  with  any  numerical  integrator  to  ev£iluate  the  period  map  and  the 
Jacobian.  Computations  for  the  DufSng  example  described  in  this  paper 
were  carried  out  using  a  Bulirsch-Stoer  integrator  adapted  from  [16]. 

Having  developed  procedures  for  computing  maps  and  Jacobians,  the 
Interpreter  can  use  Newton-Raphson  iteration  to  search  for  periodic  points 
of  a  map  /.  Naunely,  given  an  approximation  i  to  a  periodic  point  of  order 
n  for  /,  attempt  to  correct  i  by  subtracting  from  it  a  small  vector  6x  such 
that  X  —  6x  will  be  equal  to  f^”\x  —  Sx).  Linearizing  this  equation  near  x 
gives  an  approximation  for  Sx  in  terms  of  the  Jacobian  Z)/^"^(x): 

/^"^(x  —  ^x)  =  X  —  Sx 
f^^\x)-Df^^\x)Sx  »  x-^x 

Sx  w  (/ -  D/(")(x))-‘ (x  - /(’‘^(x)) 

In  terms  of  procedures  that  compute  /  and  Df,  this  computation  can  be 
accomplished  by  computing  the  orbit  of  x: 

Xo  =  X,  Xi  =  /(x,_i)  i  =  1  to  n  (2) 

and  using  the  recurrence 

Z7/I'>(x)  =  D/(xi.,)Df«-'\x) 


dSn 

dt 

dS\2 

dt 

<f^21 

dt 

dS22 

~dr 


(3) 
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(dalinc-natvork  duffing 
((x2  ■tnt«-v«Lriabl«) 

(xl  >tat«-varl&bl«) 

(p  paranatar) 

(k  paraaatar)} 

(constraint  ‘ (>  ,(rata  xl)  .x2)) 

(constraint  ‘ («  , (rata  x2) 

(-  (*  ,p  (cos  t))  (+  (*  ,k  ,x2)  (*  ,xl  ,xl  ,xl)))))) 


(lambda  (k  p) 

(lambda  (astata*) 

(lat  ((t  (vactor-raf  astataa  0)) 

(xl  (vactor-raf  astata*  1)) 

(x2  (vactor-raf  aatataa  2))) 

(vactor  1 
x2 

(♦  (a  -1  k  x2)  (a  -1  xl  xl  xl)  (a  p  (cos  t))))))) 

(lambda  (k  p) 

(lat  ((gl  (a  -1  k))) 

(lambda  (avarstataa) 

(lat  ((t  (vactor-raf  avarstataa  0)) 

(xl  (vactor-raf  avarstataa  l)) 

(x2  (vactor-raf  avarstataa  2)) 

(dal.xl.xl  (vactor-raf  avarstataa  3)) 

(dal.xl.x2  (vactor-raf  avarstataa  4)) 

(dal.x2.xl  (vactor-raf  avarstataa  8)) 

(dal.x2.x2  (vactor-raf  avarstataa  6))) 


(lat  ((g2  (a  xl  xl))) 
(lat  ((g3  (a  -3  g2))) 


(vactor  1 
x2 

(+  (a  -1  g2  xl)  (a  gl  x2)  (a  p  (cos  t))) 


dal.xl.x2 

(♦  (a  gl  dal.xl.x2)  (a  g3  dal.xl.xl)) 
dal.x2.x2 

(+  (a  gl  dal.x2.x2)  (a  g3  dal.x2.xl))))))))) 


Figure  9:  The  Interpreter  manipulates  the  symbolic  expressions  that  describe  the  system, 
in  order  to  generate  numerical  procedures  for  locating  periodic  points.  Shown  here  is 
the  actual  input  that  describes  the  OufBng  system  of  equation  (1).  Below  this  is  a  Lisp 
procedure  that  evolves  the  system  state,  which  the  Interpreter  has  generated  and  compiled; 
and  below  this  is  the  procedure  compiled  to  evolve  the  associated  variational  system. 
Notice  that  the  compiler  has  performed  some  common-subexpression  removal  here  (gl, 
g2,  and  g3)  to  make  the  computation  more  efficient. 


21 


to  compute  Df^''\x). 

Given  an  initial  approximation  i  to  a  periodic  point  one  repeatedly  up¬ 
dates  X  to  X  —  Sx,  looking  for  convergence.  Since  Newton’s  method  normally 
converges  rapidly  if  it  converges  at  all,  and  the  update  procedure  is  expen¬ 
sive  (computing  Jacobians  can  require  integration)  the  Interpreter  is  quick  to 
abandon  an  attempt  to  find  a  periodic  point  if  the  sequence  does  not  appear 
to  be  converging.  In  particular,  if  the  second  correction  Sx^"^^  is  larger  than 
the  first  correction  6x,  the  Interpreter  will  conclude  that  the  initial  guess 
for  the  periodic  point  was  bad,  and  stop  the  search.^^  Newton’s  method 
will  also  fail  if  /  —  Df^”'\x)  becomes  singular.  In  addition,  the  Interpreter’s 
Newton’s-method  procedure  can  be  called  with  an  arbitrwy  predicate,  which 
successive  guesses  must  satisfy  in  order  for  the  search  to  continue.  For  exam¬ 
ple,  guesses  must  normally  remain  inside  the  bounded  region  in  state  space 
that  the  Interpreter  is  exploring.  Also,  some  ot  tne  Interpreter’s  bifurcation 
identification  procedures  require  finding  periodic  points  that  lie  within  cer¬ 
tain  specified  neighborhoods;  these  procedures  call  Newton’s  method  with 
appropriate  predicates  to  enforce  these  constraints. 

Finally,  if  the  the  sequence  of  approximations  converges  to  a  periodic 
point  X  of  the  desired  order  n,  the  Interpreter  checks  the  orbit  of  x  to  see 
whether  x  in  fact  has  order  smaller  than  n. 

Sensitivity  derivatives 

In  addition  to  approximating  periodic  points  in  a  parametrized  feunily  of 
maps,  the  Interpreter  also  tracks  the  loci  of  the  periodic  points  as  the  pa¬ 
rameter  varies.  This  requires  computing  the  infinitesimal  change  in  position 
of  a  periodic  point  x  of  /  with  respect  to  p.  To  see  how  to  compute  this, 
consider  first  the  case  where  x  is  a  fixed  point  of  the  map  x  /(p,  x).  For 
an  incremental  change  Sp,  compute  the  corresponding  incremental  change  Sx 
such  that  X  6x  approximates  /(p  +  Sp,x  +  Sx)  to  first  order: 

X  +  Sx  =  f(p  +  Sp,x  +  Sx) 

»  f{p,  x)  +  D^fip,  x)Sx  -f  Dpf{p,  x)Sp 
=  X  +  Dfj,{x)Sx  +  Dpfp{x)Sp 

where  Dfp(x)  is  Jacobian  at  i  of  the  map  fp-.  x  f{p,  x)  and  and  Dpfp(x) 
is  the  derivative  with  respect  to  p  of  the  (vector- valued)  function  p  /(p,  x). 

^^Section  3.2  explains  how  the  Interpreter  selects  initial  approximations. 

'^his  technique  was  suggested  to  me  by  Matthew  Halfant. 
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Subtracting  i  from  both  sides  of  this  equation  and  solving  for  Sx  gives 

«x  K  (/  -  D/,(i))-'DMx)6p  (4) 

If  fp  is  given  by  explicit  equatiuos  (the  case  of  a  discrete  dynamical  sys 
tern),  the  Bifurcation  Interpreter  can  symbolically  differentiate  the  expres¬ 
sions  for  fp  with  respect  to  p  and  construct  a  numerical  procedure  that  com¬ 
putes  Dpfp{x).  For  a  system  of  ordinary  differential  equations  i  =  Hp{x,t), 
where  fp  is  the  period  map,  Dpfp{x)  can  be  computed  by  integrating  the 
components  of  dHfdp  over  the  trajectory  starting  from  x.  In  the  case  of  the 
DuflBng  system  (1)  for  example,  the  components  of  dH/dp  to  be  integrated 
are  obtained  by  differentiating  (1)  with  respect  to  p  to  obtain 


d  Uxi 
dt\dp 


dt\  dp  ) 


cos  I  —  K— - 3*;— r- 

dp  dp 


Just  as  in  computing  Jacobians,  the  Interpreter  uses  symbolic  algebra  to 
perform  this  differentiation.  It  automatically  compiles  a  sensitivity-derivative 
which  can  be  integrated  to  compute  Dpfp{x)}* 

Tracking  an  order-n  periodic  point  of  /  can  be  treated  as  tracking  a  fixed 
point  of  /("I.  However,  as  with  locating  fixed  points,  it  is  convenient  to  use 
the  recurrence  scheme  (2),  (3)  and  to  compute  Dpf^'^^x)  recursively  in  terms 
of  Dpfp{x)  by  means  of  the  relation 

This  equation  follows  from  the  chain  rule  for  partial  derivatives 
Dp{hog){x)  =  Dph{g{x))  +  Dh{g{x))Dpg{x) 
taking  h  as  fp  and  g  as 

In  summary,  whenever  the  Interpreter  successfully  applies  Newton’s  method 
to  produce  a  periodic  point  i  of  order  n,  it  saves  x,  together  with  its  orbit,  the 
value  of  p,  the  Jacobiam  Df^\x)y  and  the  sensitivity  derivative  Dpf^\x). 


‘^This  technique  for  computing  sensitivities  of  fixed  points  was  also  used  in  [1]  for 
analyzing  electrical  networks. 
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3.2  Choosing  an  initial  set  of  periodic  points 

After  compiling  the  numerical  procedures  aa  described  above,  the  Interpreter 
computes  an  initial  set  of  periodic  points  that  it  will  later  extend  to  families. 
The  parameter  interval  to  be  explored  is  divided  into  equal-length  subinter¬ 
vals  at  points  p,.  For  each  p,,  the  Interpreter  attempts  to  find  the  stable 
periodic  points  of  /p, .  The  algorithm  for  finding  these  periodic  points  pro¬ 
duces,  at  each  p,  ,  a  collection  of  periodic  points  «ind  possibly  some  “unknown 
orbits,”  which  cirise  from  trajectories  that  do  not  lead  to  periodic  points.^® 
The  Interpreter  next  compares  the  collection  of  points  discovered  at  adjzw:ent 
values  Pi,  Pi^\.  If  these  are  not  equivadent — that  is,  if  the  two  collections  do 
not  contain  the  same  number  of  periodic  points  with  corresponding  orders, 
and  the  same  number  of  unknown  orbits — then  the  interval  from  p,-  to  pi^i  is 
bisected,  and  periodic  points  are  computed  at  the  midpoint.  Bisection  con¬ 
tinues  until  adjacent  values  of  p  either  have  equivalent  periodic-point  sets,  or 
else  are  closer  than  some  pre-specified  minimum  distance. 

To  compute  periodic  points  at  a  parameter  vadue  p,  the  Interpreter  first 
guesses  approximations  for  the  periodic  points  by  means  of  the  “unravelling 
algorithm”  of  Hsu  and  Guttalu  (8,  9).  This  begins  by  imposing  a  grid  on  the 
region  in  state-space  to  be  examined.  Define  a  cell-to-ceU  transformation  by 
mapping  each  grid  cell  C  to  the  cell  containing  the  point  obtained  by  applying 
/p  to  the  midpoint  of  C,  mapping  C  to  a  distinguished  value  “infinity”  if  the 
image  of  the  midpoint  under  /p  lies  the  grid.  Regarding  the  ceUs  as  the  nodes 
of  a  graph,  where  each  cell  is  linked  to  its  image  under  the  cellular  map, 
produces  a  directed  graph  in  which  each  node  has  at  most  one  successor. 
Use  a  marking  algorithm  to  traverse  the  graph  and  pau’tition  the  cells  into 
connected  components  of  the  graph.  During  the  marking,  whenever  you 
encounter  a  cycle  of  order  n  starting  at  a  cell  C,  choose  the  midpoint  of  C 
as  a  guess  for  an  order-n  periodic  point  of  /p.'® 

'®Th«*e  unknown  orbits  may  correspond,  in  actuality,  to  orbits  that  are  indeed  non¬ 
periodic,  or  to  orbits  that  have  extremely  large  period,  or  possibly  to  bad  numerical 
properties  of  /p.  For  example,  iterating  fp  for  p  extremely  cloee  to  a  bifurcation  value  may 
lead  to  unknown  orbits. 

**Note  that  any  path  through  the  graph  will  end  in  a  cycle,  a  fixed  point  (cycle  of  order 
1)  or  else  infinity.  The  connected  components  of  the  graph  correspond  to  the  “basins  of 
attraction”  of  these  cycles.  Regarding  the  corresponding  sets  of  cells  as  approximations  to 
the  actual  basins  of  attraction  of  /  leads  to  interesting  computational  methods  for  explor¬ 
ing  dynamical  systems.  The  are  described  in  [8],  which  also  presents  more  sophisticated 
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The  Interpreter  now  tries  <iach  guess  as  the  starting  point  for  a  Newtou- 
Raphson  iteration.  If  the  iteration  succeeds,  the  periodic  point  and  its  orbit 
axe  recorded.  If  Newton’s  method  fails  (e.g.,  the  successive  iterates  may 
become  unbounded,  or  the  resulting  limit  point  may  be  unstable,  or  the 
Jacobian  may  become  singular)  the  Interpreter  uses  the  guess  as  a  starting 
point  for  a  “trajectory  search,”  which  proceeds  as  follows: 

Given  an  initial  point  x  =  xq,  generate  the  trajectory  ii  =  /p(io),  X2  = 
/p(ii),  ....  As  each  new  point  Xa  is  generated,  check  to  see  if  it  is  close 
to  a  point  in  a  previously-generated  trajectory,  in  which  case  abandon  the 
search  starting  at  x.  Alternatively,  if  Xa  is  dose  to  a  point  in  the  current 
trajectory,  scan  backward  along  the  trajectory  to  find  the  smallest  n  such 
that  Xa-n  is  close  to  x^-^^  Now  generate  a  few  more  points  starting  from  x^, 
and  check  whether  the  points  lo+i,  1^+2^  remadn  close,  respectively,  to 
Xa+i_„,  Xa+2-n,  •  ■  •  •  If  SO,  usc  x*  as  the  starting  point  for  a  Newton-Raphson 
search  for  a  stable  periodic  point  of  order  n.  If  the  iteration  does  not  succeed, 
or  if  Xa  is  not  close  to  any  previous  trajectory  point,  then  continue  evolving 
the  trajectory.  If  the  number  of  points  in  a  trajectory  grows  to  exceed  some 
majcimum,  declare  the  trajectory  to  be  of  type  “unknown”  and  abandon  it. 

As  an  example,  in  producing  the  report  on  the  Duffing  system  in  sec¬ 
tion  1.1,  the  Interpreter’s  parameters  were  set  so  that  it  initially  divided  the 
parameter  interval  from  1  to  25  into  subintervals  of  size  0.5.  The  minimum 
size  for  an  interval  (at  which  bisection  ceased)  was  0.1.  The  state-space  re¬ 
gion  — 5  <  X  <  5,-10  <  j/  <  10  was  partitioned  by  a  grid  of  size  0.5  on 
a  side.  With  these  settings,  the  Interpreter  ended  up  scanning  for  periodic 
points  at  93  values  of  p.  In  all,  the  Hsu-Guttalu  algorithm  produced  307 
guesses  for  periodic  points.  Of  these,  114  were  successfully  refined  to  peri¬ 
odic  points  by  Newton-Raphson;  the  other  193  became  starting  points  for 
trajectory  se^Lrches.  At  many  values  of  p  several  initial  guesses  were  refined 
to  the  same  periodic  point,  and  after  removing  duplicates,  the  Interpreter 
had  in  all  identified  117  distinct  periodic  points  and  15  unknown  orbits.** 

variations  on  the  unravelling  algorithni  described  above. 

*^TVajectory  points  are  stored  in  a  quad  tree  [15]  and  also  linked  in  a  doubly-linked  list 
to  facilitate  finding  close  pairs  of  points  and  scanning  forwards  and  backwards  through 
the  trajectory. 

'*Note  that  the  grid  size  used  here  was  rather  large,  so  one  should  expect  the  Hsu- 
Guttalu  algorithm  to  produce  many  bad  guesses  for  periodic  points.  Choosing  a  finer 
grid,  however,  would  be  computationally  more  expensive.  More  work  needs  to  be  done  to 
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3.3  Extending  periodic  points  to  families 

Each  stable  periodic  point  now  becomes  the  “seed  point”  for  a  family  of 
periodic  points  whose  position  varies  with  p.  The  Interpreter  tracks  the 
point,  extending  each  family  as  fax  as  possible  in  the  directions  of  increasing 
zmd  decreasing  p.  During  the  tracking  process,  families  may  merge,  or  they 
may  terminate  in  bifurcations,  which  the  Interpreter  attempts  to  classify. 

To  trzw:k  a  stable  periodic  point  x  from  p,  to  pt,  the  Interpreter  chooses 
a  smaU  increment  6p  (respectively,  a  small  decrement  6p  if  p,  >  pt)  and 
computes  the  corresponding  first-order  change  8x  using  equation  (4),  together 
with  the  values  of  the  Jacobian  and  the  sensitivity  derivative  that  were  stored 
along  with  x,  as  explained  in  section  3.1.  The  Interpreter  then  uses  Newton- 
Raphson  iteration  to  attempt  to  correct  z  -I-  to  a  stable  periodic  point  of 
/  at  p,  -I-  Sp. 

In  tracking  a  periodic  point  at  discrete  steps  like  this,  one  must  be  careful 
not  to  mistakenly  “jump”  to  the  path  of  a  neau'by  periodic  point.  The  Inter¬ 
preter  attempts  to  guard  against  this  by  adaptively  adjusting  the  stepsize  Sp, 
always  choosing  it  small  enough  so  that:  (1)  the  magnitude  of  the  computed 
6x  is  small;  (2)  the  actual  periodic  point  discovered  by  Newton’s  method  is 
close  to  the  first-order  prediction  x  -f  Sx-,  and  (3)  the  change  in  magnitude  of 
the  eigenvalues  of  the  Jacobian  from  i  to  the  new  periodic  point  is  small. 

each  step,  Newton’s  method  may  fail  to  produce  a  stable  periodic 
point.  The  iteration  may  not  converge,  or  the  matrix  /  —  Dpf^''\x)  may 
become  singular,  so  that  no  periodic  point  will  be  produced.  Alternatively, 
the  periodic  point  found  may  be  unstable  (which  can  be  determined  by  ex¬ 
amining  the  eigenvalues  of  the  Jacobian),  Or,  the  point  may  be  stable,  but 
with  an  order  less  than  the  order  of  x.  Upon  such  failures,  the  Interpreter 
decreases  the  size  of  Sp  and  tries  again.  When  Sp  becomes  less  than  some 
prespecified  lower  bound,  and  Newton’s  method  still  fails,  the  Interpreter 
regards  the  failure  as  evidence  of  a  bifurcation  at  a  piirameter  value  between 
p,  and  p,  -f-  Sp.  It  records  the  nature  of  the  faiilure  and  uses  this  information 
in  am  attempt  to  claissify  the  bifurcation,  ais  described  below  in  section  3.4. 
If  Newton-Raphson  succeeds,  the  Interpreter  increaises  the  stepsize  Sp  and 
continues  tracking  towaxd  p^. 

investigate  the  tradeofb  here,  ud  to  find  good  ways  to  pick  the  initial  grid  size.  Siniilarly, 
it  is  not  obvious  how  to  choose  good  sizes  for  the  parameter  intervals.  Choosing  intervals 
that  are  large  risks  overlooking  significant  families  of  periodic  points  and  bifurcations. 


26 


The  Interprjher  cn'jpioys  this  tracking  method  to  extend  each  family  as  far 
as  possible — until  I  he  family  terminates  in  a  bifurcation,  or  the  periodic  point 
leaves  the  bounded  region  in  state-space  that  the  Interpreter  was  directed  to 
explore,  or  the  tracker  reaches  an  endpoint  of  the  parameter  interval.  As  each 
faunily  is  generated,  it  is  represented  as  a  list  of  periodic  points  in  order  of 
increasing  p.  Let  the  current  starting  value  of  a  family  be  the  smallest  vaJue  of 
p  to  which  the  family  has  so  far  been  extended,  and  the  first  knovm  instance 
is  the  periodic  point  at  that  value  of  p.  Similarly,  the  current  ending  value 
and  the  last  known  instance  give  the  parameter  and  the  periodic  point  at  the 
largest  value  of  p  to  which  the  family  has  so  far  been  extended.  A  family  is 
said  to  have  its  vanishing  known  (respectively,  its  appearance  known)  if  it  haa 
been  extended  as  far  as  possible  in  the  direction  of  increasing  p  (respectively, 
decreasing  p). 

While  tracking,  the  Interpreter  must  also  check  for  merges  among  faimilies. 
This  requires  some  care,  because  individual  families  are  tracked  sepairately, 
and  the  points  of  different  families  will  in  general  be  computed  at  different 
parameter  values.  Specifically,  the  Interpreter  proceeds  as  follows:^* 

The  parameter  interval  is  pairtitioned  at  a  number  of  equad-spaced  syn¬ 
chronization  values.  The  Interpreter  chooses  a  family,  say  family  A,  whose 
varnishing  is  not  yet  known.  Let  x  be  the  first  known  instance  of  A,  at 
paurameter  value  p,.  Choose  pt  to  be  the  minimum  of; 


1.  the  right  most  endpoint  of  the  parameter  interval 

2.  the  smallest  value  after  p,  that  is  the  current  starting  value  of  a  family 
whose  appearance  is  not  yet  known 

3.  the  smaUest  synchronization  value  lairger  than  p, 

4.  the  smallest  value  larger  than  p,  at  which  there  is  a  bifurcation 


If  the  tracker  fails  to  extend  A  as  far  as  pt,  then  the  type  of  failure  (point 
going  out  of  bounds  or  reaching  a  bifurcation)  shows  how  A  vanishes,  amd  the 
Interpreter  accordingly  maurks  the  A’s  vanishing  ais  known.  If  the  extension 
to  Pt  is  successful,  then  the  the  Interpreter  checks  for  merges  of  A  with  other 
families.  There  are  several  cases  to  consider,  depending  on  how  p*  wais  chosen: 

'®For  simplicity,  we  describe  the  tracking  process  in  the  direction  of  increasing  p.  Track¬ 
ing  in  the  direction  of  decreasing  p  proceeds  analogously. 
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Figure  10:  Merging  families:  (a)  If  family  A  extends  to  meet  the  first  known  Lnstance 
of  B,  the  two  families  are  merged  to  form  a  single  family  AB.  (b)  If  A  coincides  with 
a  point  in  the  middle  of  B  (at  a  synchronization  value  of  the  parameter),  B  is  split  into 
two  families  B\  and  B^.  A  and  Bx  vanish  together  at  a  bifurcation  somewhere  before 
the  synchronization  point,  and  B^  appears  there.  Note  that,  in  this  figure,  the  horizontal 
direction  represents  the  p-axis,  while  the  vertical  direction  represents  the  location  the 
points  in  state  space  (projected  onto  a  one-dimensional  subset). 

(1)  li  pi  is  the  endpoint  of  the  parameter  interval,  then  the  Interpreter 
marks  A  as  persisting  to  the  end  of  the  interval,  and  the  extension  of  A  is 
complete. 

(2)  If  pt  is  the  current  starting  value  of  a  family  B,  the  Interpreter  com¬ 
pares  the  currently  tracked  point  with  the  first  known  instance  of  B.  If  the 
points  are  the  same,  then  A  and  B  are  merged  as  shown  in  figure  10a.  If  not, 
tracking  of  A  continues  from  p*. 

(3)  If  Pt  is  a  synchronization  value,  the  currently  treu:ked  point  is  compared 
with  the  points  of  other  families  pt.  If  it  matches  a  point  of  some  other  family 
B,  then  A  amd  B  must  merge  at  some  value  between  p,  and  pi.  This  merge 
corresponds  to  a  bifurcation;  and  B  must  actuadly  be  two  families:  Bi  (the 
part  of  B  before  the  bifurcation)  amd  B^  (the  part  of  B  after  the  bifurcation). 
At  the  bifurcation,  A  and  Bi  vanish  and  appears.  The  Interpreter  searches 
the  interval  between  p,  and  pt  to  locate  the  merge  point  and  replaces  B  by 
Bi  and  B-x  (see  figure  10b). 

Locating  the  precise  merge  point  here  is  slightly  tricky.  Families  A  and 
B  coincide  at  p*  and  are  distinct  at  p,,  so  the  merge  point  can  be  found 
by  binary  search.  However,  the  periodic  points  on  A  and  B  have  not  been 
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computed  at  the  same  parameter  values  between  p,  2uid  pt  and  so  cannot  be 
compared  directly.  Comparing  the  families  at  intermediate  parameter  values 
thus  requires  computing  additional  points  ol  A  or  B  (or  both).  This  is 
accomplished  by  tracking  from  existing  points  at  nearby  parauneter  values.*® 

Finally  (4),  if  pt  is  a  bifurcation  value,  the  Interpreter  checks  to  see  if  A 
vanishes  at  the  bifurcation.  If  not,  tracking  of  A  continues  at  P(.  Otherwise, 
A  is  compared  with  the  other  families  currently  known  to  vanish  at  the 
bifurcation.  (This  comparison  ma>  require  generating  additional  periodic 
points,  as  in  the  previous  paragraph.)  If  A  does  coincide  with  soms  family 
B  that  has  already  been  tracked  to  the  bifurcation  at  pt,  then  A  must  in  fact 
merge  with  B  at  another  bifurcation  before  pt,  as  in  case  (3).  Alternatively, 
if  i4  is  a  new  family  vanishing  at  the  bifurcation,  then  the  Interpreter  marks 
it  as  such. 

In  this  way,  tracking  proc.eed3  until  all  faunilies  have  been  extended  so 
that  their  appearance  and  vanishing  is  known. 

3.4  Identifying  bifurcations 

When  the  periodic-point  tracking  algorithm  fails,  the  Interpreter  assumes 
that  it  has  encountered  a  bifurcation,  and  it  attempts  to  classify  the  bifur¬ 
cation.  There  are  three  ways  in  which  tracking  a  stable  periodic  point  might 
fail:  (1)  Newton-Raphson  iteration  fails  to  converge,  so  that  no  new  periodic 
point  is  found  at  the  next  parameter  value;  (2)  the  iteration  converges,  but 
to  an  unstable  periodic  point;  (3)  the  iteration  converges  to  a  stable  periodic 
point,  but  of  lower  order  than  the  point  being  tracked. 

Comparing  these  failure  types  against  the  descriptions  of  the  generic  bi¬ 
furcations  given  in  section  2.1  and  figure  7  indicates  how  to  begin  matching 
tracking  failures  ageunst  bifurcation  types.  In  case  (3),  for  example,  the  only 
bifurcation  in  which  a  stable  periodic  point  transforms  into  a  stable  periodic 
point  of  a  lower  order  is  the  supercritical  flip,  where  a  point  of  order  2n 
changes  to  a  point  of  order  n.*’  Accordingly,  if  the  Interpreter  encounters 

*®There  is  also  the  possibility  that  this  secondary  tracking  may  fail,  since,  after  all,  one 
is  close  to  a  bifurcation,  this  case  the  Interpreter  assumes  that  a  tracking  failure  here  is 
due  to  the  bifurcation  it  is  searching  for. 

**  Notice,  with  reference  to  figure  7,  that  this  transition  corresponds  to  crossing  through 
the  supercritical  flip  bifurcation  from  right  to  left.  The  Interpreter  must  recognise  bifur¬ 
cations  encountered  from  either  direction. 
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a  failure  of  type  (3),  it  expects  that  the  order  of  the  new  point  will  be  half 
the  order  of  the  old  point,  and  that  an  eigenvalue  of  map  is  crossing  the 
unit  circle  at  —1.  If  both  these  conditions  hold,  the  Interpreter  will  attempt 
to  identify  the  transition  as  a  supercritical  flip.  Otherwise,  it  will  report  a 
bifurcation  of  unknown  type. 

For  the  Interpreter,  “identifying  a  bifurcation”  means  finding  the  com¬ 
plete  complement  of  stable  and  unstable  periodic  points  that  should  coalesce 
at  the  bifurcation.  Thus,  at  a  supercritic^d  flip,  there  is  a  stable  point  of 
order  n  to  one  side  of  the  bifurcation,  and,  to  the  other  side,  a  stable  point 
of  order  2n  and  an  unstable  point  of  order  n.  If  the  Interpreter  has  found  the 
two  stable  points,  it  attempts  to  locate  the  expected  unstable  point,  using 
methods  that  are  descirbed  below.  If  it  succeeds,  it  reports  the  bifurcation 
as  a  supercritical  flip.  Otherwise  it  reports  the  bifurcation  ais  unknown. 

In  general,  the  Interpreter  begins  its  attempt  to  claissify  a  bifurcation 
by  examining  the  eigenvalues  of  the  point  near  the  tracking  failure.  If  an 
eigenvalue  crosses  the  unit  circle  near  —1,  the  bifurcation  (if  it  is  one  of  the 
recognized  types)  must  be  a  supercritical  or  a  subcritical  flip;  crossing  the 
unit  circle  near  1  leads  to  a  fold,  a  transcritical  bifurcation,  a  supercritical 
pitchfork,  or  a  subcritical  pitchfork;  crossing  the  unit  circle  off  the  real  axis 
leads  to  a  supercritical  or  subcritical  Niemark  bifurcation.  In  each  case,  the 
Interpreter  scfirches  near  the  bifurcation  for  additional  stable  and  unstable 
periodic  points,  until  it  flnds  enough  points  to  match  the  bifurcation  against 
one  of  the  known  types.^* 

If  a  recognized  local  portrdt  is  found — i.e.,  if  the  number  of  stable  and 
unstable  periodic  points  at  the  bifurcation  matches  the  pattern  for  one  of  the 
known  bifurcation  types — the  Interpreter  records  a  successful  local  identifi¬ 
cation.  Otherwise,  the  Interpreter  records  the  bifurcation  type  as  unknown. 
In  either  case,  the  result  of  the  cljissification,  together  with  the  set  of  periodic 
points  found  by  the  search  is  passed  to  the  next  phase  of  the  program,  which 
criticizes  these  local  identifications  in  the  light  of  more  global  information, 
as  discussed  in  section  3.5. 

the  case  of  a  Niemark  bifurcation,  there  is  a  stable  fixed  point  to  one  side  of  the 
bifurcation,  and,  to  the  other  side,  an  unstable  fixed  point  and  a  limit  cycle — a  stable  limit 
cycle  for  a  supercritical  Niemark  and  an  unstable  limit  cycle  for  a  subcritical  Niemark. 
Currently,  the  Interpreter  looks  only  for  fixed  points,  not  limit  cycles.  Thus  it  cannot 
distinguish  a  supercritical  Niemark  from  a  subcritical  Niemark,  and  reports  any  such 
bifrircation  simply  as  a  “Niemark  bifurcation.” 
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Finally,  any  nev/  stable  periodic  points  located  during  this  process  become 
seed  points  for  new  famihec,  which  must  be  tracked  in  the  direction  away  from 
the  bifurcation.  These  are  added  to  the  collection  of  families  to  be  extended 
using  the  method  of  section  3.3,  and  the  process  continues. 

Searching  for  periodic  points  near  bifurcations 

When  the  Interpreter  searches  for  a  new  periodic  point  near  a  bifurcation,  it 
guesses  a  location  for  the  point,  and  uses  this  guess  as  the  starting  point  for 
a  Newton- Raphson  iteration,  attempting  to  find  the  actual  new  point. 

The  guess  for  the  new  periodic  point  depends  upon  the  configuration  of 
periodic  points  already  observed  near  the  bifurcation.  As  figure  7  indicates, 
when  one  is  very  close  to  the  bifurcation,  the  periodic  points  tend  to  approach 
it  along  either  parabolic  or  straight-line  paths.*^  The  Interpreter  uses  this 
local  geometry  to  extrapolate  locations  for  the  unknown  periodic  points. 

For  example,  suppose  that  the  Interpreter  is  investigating  a  bifurcation  at 
which  an  eigenvalue  is  approaching  1,  and  that  there  are  two  stable  periodic 
points  (of  the  same  order)  to  one  side  of  the  bifurcation  and  one  stable 
periodic  point  to  the  other  side.  Given  the  possibilities  in  figure  7,  this  ought 
to  be  a  supercritical  pitchfork  bifurcation,  and  there  ought  to  be  another 
unstable  periodic  point  on  the  same  side  as  the  two  stable  points.  Moreover, 
the  local  geometry  near  a  supercritical  pitchfork  indicates  that  the  two  stable 
points  are  located  approximately  symmetrically  with  respect  to  the  unstable 
point.  The  Interpreter  therefore  searches  for  an  unstable  periodic  point  at 
the  same  parameter  value  as  the  two  stable  points,  and  at  their  average 
position  in  state  space.  This  particular  search  method  is  called  averaging. 
The  same  method  is  used  when  there  are  two  unstable  points  to  one  side 
of  the  bifurcation  and  an  stable  point  to  the  other  side  to  search  for  a  new 
stable  point  at  the  average  position  of  the  unstable  points. 

In  general,  each  of  the  Interpreter’s  search  methods  is  triggered  by  a 
particular  pattern  of  stable  and  unstable  periodic  points.  To  investigate  a 
bifurcation,  the  Interpreter  tries  each  applicable  search  method,  looking  for 
new  periodic  points.  As  new  periodic  points  are  discovered,  different  search 
methods  become  applicable.  If  all  applicable  search  methods  have  been  tried, 

^^These  local  geometries  can  be  verified  by  truncating  the  power-series  expansion  for 
f{p,x)  near  the  bifurcation  and  solving  for  the  paths  of  the  periodic  points.  This  is  in 
fact  the  basic  technique  for  cataloguing  the  types  of  generic  bifurcations.  See  [19]  or  [6] 
for  details. 
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and  the  pattern  of  fixed  points  is  not  recognized,  the  Interpreter  reports  the 
bifurcation  as  unknown. 

Here  is  the  rest  of  the  search  methods  used  at  bifurcations  where  an 
eigenvalue  approaches  1.  For  these  bifurcations  the  periodic  points  all  have 
the  same  order.  These  methods  are  also  illustrated  schematically  in  figure  11. 

•  re/Iect. Pattern — a  stable  and  an  unstable  periodic  point  to  one  side  of 
the  bifurcation.  Search  for  a  second  stable  point  at  the  same  parameter 
value  by  reflecting  the  stable  point  in  the  unstable  point.  (Analogously, 
search  for  a  second  unstable  point  by  reflecting  the  unstable  point  in 
the  stable  point.) 

•  extrapolate  across:  Pattern — a  stable  point  and  an  unstable  point  to 
one  side  of  the  bifurcation,  say,  at  parameter  value  po,  and  an  unstable 
point  to  the  other  side  at  pb.  Search  for  a  stable  point  at  pb  as  follows: 
At  Pa  And  the  vector  from  the  unstable  point  to  the  stable  point.  At 
Pb  subtract  this  vector  from  the  unstable  point,  and  search  for  a  stable 
point  here.  Analogously,  given  instead  a  stable  point  at  pb,  search  for 
an  unstable  point  by  extrapolating  from  the  stable  and  unstable  points 
at  Pa. 

•  extrapolate  through:  Pattern — two  stable  points  and  one  unstable  point 
to  one  side  of  the  bifurcation.  Search  for  a  stable  point  on  the  other 
side,  at  the  same  position  in  state  space  as  the  unstable  point. 

•  quadratic  extrapolate:  Pattern — a  single  stable  point  to  one  side  of  the 
bifurcation  at  pb.  Search  for  an  unstable  point  also  at  pb,  using  the 
fact  that,  for  folds  or  pitchforks,  points  very  close  to  the  bifurcation 
approach  it  along  parabolic  paths.  The  method  requires  knowing  not 
only  the  location  of  the  point  at  pb,  but  adso  the  periodic-point  found 
by  the  tracker  at  the  previous  step  to  pb,  say  at  parameter  vaJue  pc- 
Let  Pa  be  the  parameter  value  where  the  tracking  failed.  Fit  a  branch 
of  a  parabola  through  the  points  (xp^,Pc),  (xpuPfc),  and  (0,pa).  Tzdce, 
as  a  guess  for  the  unstable  point,  the  point  on  the  opposite  branch  of 
the  parabola  at  p(,. 
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Figure  1 1 :  These  figures  indicate  the  methods  used  by  the  Interpreter  to  search  for  new 
stable  and  unstable  periodic  points  near  a  bifurcation  where  an  eigenvalue  approaches  1. 
Solid  dots  indicate  stable  points  and  hollow  dots  indicate  unstable  points.  The  coordinate 
system  drawn  here  has  the  bifurcation  at  the  origin.  The  horizontal  axis  gives  the  param¬ 
eter  direction  p.  The  vertical  axis  represents  the  two-dimensional  state  space.  The  curves 
shown  for  the  extrapolation  methods  indicate  the  assumed  linear  or  parabolic  paths  of  the 
peiiudic  pouits  near  liie  bifurcation,  which  form  the  basis  for  the  extrapolation.  Compare 
this  figure  with  the  bifurcation  diagrams  in  figure  7  to  see  the  full  local  geometries  that 
these  search  methods  are  attempting  to  complete. 
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Here  is  an  example,  taken  from  a  test  of  the  Interpreter,  which  illustrates 
these  search  methods  in  action.  The  test  map  here  is 

fp  ••  (xux^)  I  ~  X 

This  map  has  a  transcritical  bifurcation  at  p  =  1,  (ii,X2)  =  (0,0),  which 
the  Interpreter  successfully  identified.  In  the  test,  the  Interpreter  first  en¬ 
countered  the  bifurcation  while  tracking  a  fixed  point  for  decreasing  p.  The 
tracker  found  a  stable  fixed  point  before  the  bifurcation  at  p  =  1.0004  and  an 
unstable  fixed  point  after  the  bifurcation  at  p  =  0.999.  Over  the  transition, 
an  eigenvalue  crossed  the  unit  circle  at  1.  The  Interpreter  fi^st  used  quadratic 
extrapolation  of  the  stable  point.  This  successfully  located  an  imstable  point 
before  the  bifurcation.  The  Interpreter  next  attempted  to  find  a  second  stai- 
ble  point  before  the  bifurcation  by  reflecting  the  stable  point  in  the  unstable 
point.  This  attempt  failed.  The  Interpreter  then  extrapolated  the  unstable 
point  across  the  bifurcation  to  produce  a  stable  point  after  the  bifurcation. 
This  resulted  in  stable-unstable  pairs  both  before  and  after  the  bifurcation, 
and  the  Interpreter  announced  this  as  a  transcritical  bifurcation. 

Index  scan 

The  search  methods  described  above  are  computationally  inexpensive,  in  that 
they  propose  a  definite  guess  for  a  period-point  location  and  apply  Newton- 
Raphson;  no  search  of  state-space  is  required.  However,  these  methods  are 
not  sufficient.  In  the  cMe  of  approaching  a  supercritical  flip  bifurcation,  for 
instance,  the  Interpreter  is  tracking  a  stable  periodic  point  of  order  n  that 
suddenly  becomes  unstable  when  an  eigenvalue  crosses  the  unit  circle  at  —1. 
There  should  be  a  stable  periodic  point  of  order  2n  on  the  other  side  of 
the  bifurcation,  but  there  is  no  obvious  guess  for  the  location  of  this  point. 
The  best  one  can  do  is  to  sccirch  a  neighborhood  of  the  bifurcation  in  the 
two-dimensional  state  space. 

Fortunately,  the  two-dimensional  search  can  be  replaced  by  two  one¬ 
dimensional  searches  (which  is  computationally  much  cheaper)  using  a  method 
described  by  Hsu  [7],  based  upon  the  Poincare  index.  If  p  is  a  map  from  the 
plane  to  the  plane  and  C  is  a  closed  curve,  define  I{g,  C) — the  index  of  g 
over  C — to  be  the  number  of  rotations  over  the  path  of  the  varying  direction 
from  X  to  g{x): 

J(s,C)  =  ^  j,^JL«(x)  (5) 
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where  v(x)  is  the  vector  pointing  from  x  to  ^(x).  If  C  is  a  simple  close-i 
curve,  then  I{g,  C)  is  equal  to  the  sum,  over  the  fixed  points  x  of  ^  enclosed 
by  C,  of  the  local  indices  I{g,  x),  where  I {g,  i)  =  I{g,  C*)  for  a  small  path 
enclosing  the  fixed  point  x  (and  no  other  fixed  points  of  gf).  Moreover,  the 
local  index  of  a  nondegenerate  fixed  point  i  is  determined  by  th*  Jaeobian 
of  g  at  X'}* 


f —1  if  det(Dg(x)  —  I)  <  0 
if  det(Dg(x)  —  /)  >  0 


(6) 


The  formulas  lead  to  an  index  scan  search  method  for  finding  periodic 
points,  essentially  as  described  in  [7].  For  example,  to  find  the  stable  periodic 
point  of  order  2n  on  the  other  side  of  a  supercritical  flip  bifurcation  of  a  map 
/,  begin  with  the  known  unstable  order-n  periodic  point  x„  of  /  and  apply 
equation  (6)  with  g  =  to  compute  the  local  index  of  at  x„.  Using 
binary  search,  attempt  to  find  circles  C,  of  radius  r  and  Cr+sr  of  radius  r  +  Sr 
such  that  the  index  around  as  computed  by  (5),  is  equal  to  the  local  index 
of  g  at  Xu,  and  the  index  around  Cr+Sr  is  not.^*  There  must  be  another  fixed 
point  of  in  the  annulus  between  Cr  and  Cr+Sr-  Now  search  for  the  fixed 

point  along  the  annulus,  choosing  each  successive  point  as  the  start  for  a 
Newton-Raphsou  iteration. 

This  search  method  is  much  more  computationally  expensive  than  the 
methods  described  in  the  previous  section.  The  Interpreter  resorts  to  index 
scanning  only  when  none  of  the  other  methods  apply,  or  when  they  have 
already  been  tried  and  yet  the  bifurcation  has  still  not  been  identified. 

Note:  In  applying  equation  (5)  the  Interpreter  does  not  integrate  angles 
directly,  but  rather  uses  an  algorithm  due  to  Lelauid  [10],  which  computes 
winding  numbers  by  counting  the  varying  direction’s  signed  o^ssings  with 
tue  positive  i-axis. 


^^These  relations  are  proved  by  Hsu  in  [7],  Note  that  the  index  here  is  not  lae  same  as 
the  tnsiabtliiy  \ndez  (number  of  eigenvalues  with  magnitude  greater  than  1),  which  also 
plays  a  role  in  dynamical  systems  theory. 

^‘Occasionally,  numerical  error  will  make  it  impossible  to  start  the  search  by  locating 
a  small  circle  for  which  the  index  is  the  same  as  the  local  index,  in  which  case  the  In¬ 
terpreter  will  abandon  the  index  scan  for  this  point.  Also,  the  Interpreter  will  not  begin 
to  investigate  a  suspected  period  doubling  if  the  order  of  the  periodic  point  exceeds  some 
predefined  maximum  (order  32  for  the  examples  in  this  paper). 
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3.5  Criticizing  local  results 

At  this  point  in  the  process,  the  Interpreter  has  isolated  a  collection  of  fam¬ 
ilies  of  periodic  points.  Each  family  has  been  fully  extended,  and  the  bi¬ 
furcations  at  which  the  families  terminate  have  been  tentatively  classified. 
These  classifications  axe  based  upon  purely  local  information — analyses  of 
individual  bifurcations,  triggered  by  the  tracking  failure  of  individual  peri¬ 
odic  points.  The  local  methods  work  well  when  the  map  fp  can  be  computed 
very  accurately.  In  dozens  of  test  cases,  where  fp  was  given  by  simple  arith¬ 
metic  formulas,  local  searches  successfully  classified  all  bifurcations.  How¬ 
ever,  when  fp  cannot  be  so  accurately  computed — as  in  the  case  of  of  period 
map  that  must  be  numerically  integrated — numerical  errors  can  confuse  the 
local  methods.  Therefore,  after  all  bifurcations  and  famihes  have  been  found, 
the  Interpreter  reexamines  its  classifications  and  attempts  to  identify  and 
correct  errors. 

The  Interpreter  begins  by  removing  families  that  could  not  be  tracked, 
l.e.,  that  contain  only  a  single  periodic  point  that  could  not  be  extended. 
Presumably  these  are  flukes  due  to  numerical  error  or  to  extremely  bad  lo¬ 
cal  behavior.  Any  bifurcation  involving  such  a  fcunily  has  its  type  reset  to 
unknown.  Two  bifurcations  that  are  connected  by  such  a  f^lmily  are  merged, 
the  type  of  the  merged  bifurcation  is  set  to  unknown. 

The  Interpreter  next  tries  to  determine  whether  two  bifurcations  that 
appeared  to  be  distinct  are  in  reality  the  same  bifurcation.  Figure  12  shows 
an  example  of  this  kind  of  error.  Bifurcations  B\  and  were  encountered 
by  tracking  stable  points  in  the  direction  of  increasing  p,  and  each  traicking 
failed  when  an  eigenvalue  approached  1.  At  each  bifurcation,  the  Interpreter 
found  an  unstable  periodic  point  before  the  bifurcation,  and  so  classified  each 
bifurcation  as  a  fold.  Also,  there  is  a  nearby  bifurcation  B3  that  was  found 
by  tracking  a  stable  periodic  point  in  ♦he  direction  of  decreasing  p  until  an 
eigenvalue  approached  1.  At  B3  none  of  the  search  methods  produced  a 
new  periodic  point,  and  the  Interpreter  classified  B3  as  unknown.  In  fact, 
however,  these  three  bifurcations  are  the  same  one.  It  is  a  supercritical 
pitchfork,  smd  the  unstable  points  at  Bi  and  B3  belong  to  the  same  family 
of  unstable  points.  Due  to  numerical  error,  the  bifurcations  were  found  at 
slightly  different  locations — locations  sufficiently  different  to  be  outside  the 
predefined  tolerance  that  the  Interpreter  uses  to  check  whether  a  point  being 
tracked  has  merged  into  a  previously-discovered  bifurcation. 
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Figure  12:  Numerical  error  can  confuse  the  Interpreter’s  local  identification.  Here, 
three  “different”  bifurcations  Bj,  and  B3,  were  discovered  by  tracking  three  different 
families  of  stable  periodic  points  to  slightly  different  locations.  In  fact,  this  is  probably 
a  single  supercritical  pitchfork  bifurcation,  and  the  two  unstable  periodic  points  at  Bi 
and  flj  belong  to  the  same  family.  The  Interpreter  reexamines  its  local  classifications  and 
attempts  to  recognize  such  errors. 

One  might  argue  that  such  errors  could  be  avoided  simply  by  choosing  a 
larger  tolerance  during  tracking.  3ut  that  would  lead,  in  other  situations,  to 
merging  nearby  bifurcations  that  are  in  fact  distinct.  Instead,  the  Interpreter 
waits  until  all  bifurcations  have  been  located  and  then  checks  for  appropriate 
merges.  Suppose  that  B  was  identified  as  a  known  type,  with  all  required 
stable  and  unstable  periodic  points  discovered.  Examining  the  list  of  known 
bifurcations  in  figure  7  shows  that  B  cannot  possibly  merge  with  another 
bifurcation  to  form  another  one  of  known  type — unless  B  was  identified  as  a 
fold,  in  which  case  B  may  actually  be  a  pitchfork  or  a  transcritical  bifurcation 
at  which  the  Interpreter’s  search  methods  failed  to  discover  all  the  local 
families. 

Thus,  the  Interpreter  reexamines  only  the  folds  and  unknown  bifurcations 
to  find  pairs  that  are  suspiciously  close.  During  the  tr^w:king  phase,  each 
bifurcation’s  parameter  value  p  wais  localized  to  an  interval  of  size  6p.  The 
Interpreter  uses  equation  (4)  to  approximate  the  chemge  6x  in  periodic-point 
location  corresponding  to  a  parauneter  change  of  size  6p.  Two  bifurcations 
are  merged  if  their  parameter  intervals  overlap  and  their  periodic  points  are 
within  a  distance  of  the  magnitude  of  6x. 

Finally,  the  Interpreter  reexamines  all  bifurcations  that  are  still  classified 
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as  unknown,  and  it  opportunistically  tries  to  identify  them  based  on  very  lib¬ 
eral  criteria.  For  example,  a  stable  family  of  order  n,  tracked  to  a  bifurcation 
with  eigenvalue  approaching  —1,  is  probably  a  flip,  even  if  the  Interpreter 
waa  unable  to  find  the  point  of  order  2n.  The  Interpreter  checks  the  bifur¬ 
cation  at  the  other  end  of  the  order-n  family,  and  if  this  is  a  supercritical 
flip,  it  marks  the  unknown  bifurcation  as  probable  supercritical  flip,  on  the 
grounds  that  chains  of  supercritical  flips  (period-doubling  cascades)  are  fairly 
common.  These  opportunistic  identifications  are  flagged  as  having  been  per¬ 
formed  at  this  stage  rather  than  during  tracking,  since  they  are  presumably 
less  reliable  than  identifications  that  have  been  confirmed  by  finding  all  the 
expected  local  periodic  points. 

3.6  Generating  the  family  report 

Once  bifurcations  and  families  have  been  isolated  and  classified,  it  is  a  simple 
matter  to  produce  a  family  report.  The  Interpreter  first  scans  the  bifurca¬ 
tions,  and  annotates  the  ones  involving  symmetric  pairs  or  period- doubling 
cascades.  Any  chain  of  two  or  more  supercritical  flips  is  marked  as  a  period¬ 
doubling  cascade.  Symmetric  paurs  are  recognized  as  evolving  from  super¬ 
critical  pitchforks  whose  two  branches  lead  to  sequences  of  bifurcations 
and  B},  where,  for  each  t,  Bf  and  Bj  have  the  same  type  and  are  located  at 
the  sarnie  value  of  p. 

The  annotated  bifurcations  and  families  can  be  viewed  ais  a  graph,  where 
the  nodes  are  the  bifurcations  and  the  arcs  awe  the  families.  The  Interpreter 
separates  these  into  connected  components  of  the  graph,  which  form  the 
claisses  of  families.  Then  it  assigns  names  to  the  classes  and  to  the  faunilies 
and  bifurcations  within  each  class,  numbered  in  the  direction  of  increaising 
p.  The  final  family  report  is  a  list,  for  each  claiss,  of  the  families  or  groups  of 
families  in  the  claiss.  The  list  for  class  A  of  the  Duffing  system  analysis,  for 
instance,  was  shown  in  section  2.2. 

The  bifurcation  data  structures  also  retaun  information  that  was  gener¬ 
ated  during  the  tracking  and  recognition  phases.  Here,  for  example,  is  the 
first  bifurcation  in  class  E  of  the  Duffing  system  (figure  1).  This  is  a  super¬ 
critical  pitchfork  at  which  family  Eo  splits  into  Eifl  and  Eiy. 

(bifurcation ;bif . 1 
(paranatar-Tanga  (21.3493  21.3603)) 

(faailiaa-Taniahlng  (iCfaailjtE.O] )) 
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(J«aili*8-*ippearing  (JfLfamilyrE.  1.0]  #[la*ily  :E.  1 . 1]  )) 

(id-inf o 

((bif-id-info  (#Cid-info3  (cannot  identify  faddl*  bifurcation))) 
(bif-id-info  (#Cid-info]  (cannot  identify  saddle  bifurcation))))) 
(coanents  ( (aerged-to-incorporate  # [bifurcation:  bif.9]))) 
(classified-by-critic  true) 

(type  supercritical-pitchfork) 

(synnetric-pair  (split  to  fora  syaaetric  pair)) 

(direction  paraaieter-increasing) ) 

The  bif-id-info  slot  here  contains  information  generated  by  the  local  classi¬ 
fication  algorithm  (section  3.4),  including  pointers  to  the  stable  and  unstable 
periodic  points  discovered  by  the  local  search,  together  with  a  comment  that 
the  local  recognition  algorithm  failed — the  Interpreter  found  that  this  was 
a  saddle  bifurcation  (eigenvalue  approaching  1)  but  was  unable  to  identify 
it.  There  are  two  bif-id  entries  here  because  the  bifurcation  was  formed 
by  merging  two  bifurcations,  each  with  a  local  identification  that  failed,  ais 
described  in  section  3.5.  As  noted,  the  final  classification  of  this  bifurcations 
as  a  supercritical  flip  was  accomplished  in  the  critic  i  :ase,  after  the  merge. 

Finally,  the  lists  of  families  and  bifurcations  ar*^  aversed  by  a  simple 
text-generator,  or  a  graphics  generator,  to  produc  e  kind  of  text  and 
graphical  output  illustrated  in  section  1.1.  Not  all  information  in  the 

lists  is  reflected  in  this  output,  but  the  information  ains  as  pau't  of  the 

family  report,  available  for  further  processing. 

4  Discussion 

The  ability  to  progress  beyond  raw  numerical  data  ;  uncover  underlying 
qualitative  phenomena  is  sometimes  called  insight.  It  is  intriguing  that  a 
degree  of  this  “insight”  can  be  achieved  automatically,  with  only  a  few  sim¬ 
ple  methods.  Although  the  Bifurcation  Interpreter  combines  techniques  from 
numerical  computing,  symbolic  algebra,  and  knowledge  encoding,  it  draws 
upon  each  of  these  areas  to  only  a  very  limited  degree.  Conversely,  con¬ 
sidering  each  of  these  areas  in  turn  immediately  suggests  extensions  and 
improvements  to  the  present  implementation. 

The  most  obvious  extension  is  to  higher-dimensional  systems.  Note  that 
the  index  scan  method  for  searching  for  periodic  points  near  bifurcations 
(section  3.4)  is  the  only  algorithm  in  the  current  implementation  that  relies  on 
the  fact  that  the  state-space  is  two-dimensional.  Additional  computational 
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methods  for  “extending  a  periodic  point  through  a  bifurcation,”  which  may 
be  successful  in  higher  dimensions,  are  discussed  by  Parker  and  Chua  [13].^® 
Exploring  systems  where  the  parameter  space  has  more  than  one  dimension 
is  more  challenging;  for  one  thing,  the  generic  bifurcations  of  such  systems 
have  not  been  mathematically  classified.  But  it  should  be  possible  to  extend 
the  Interpreter  to  map  out  these  higher-dimensional  families  as  well,  amd  to 
identify  the  bifurcations  encountered,  at  least  along  codimension- 1  subspaces 
of  the  parameter  space. 

In  the  area  of  improving  the  numerical  methods,  one  straiightforward 
extension  is  to  have  the  Interpreter  track  and  report  on  families  of  unstable 
periodic  points,  in  addition  to  the  stable  ones.  This  would  produce  a  more 
complete  picture  of  the  dynamical  system.  A  second  improvement,  in  the 
case  of  ordinary  differential  equations,  is  to  have  the  prograun  use  information 
based  upon  the  complete  trajectories  of  points  rather  tham  just  the  period 
map.  For  instance,  in  testing  for  a  symmetric  pair  of  families,  the  Interpreter 
should  verify  the  symmetry  of  the  trajectories.  In  addition,  much  work  needs 
to  be  done  in  choosing  the  tolerances  of  the  numerical  routines — such  as  the 
minimum  distance  at  which  two  bifurcations  are  considered  to  be  distinct. 
One  possible  approach  here  is  to  run  the  program  repeatedly  with  different 
tolerances  until  the  high-level  qualitative  report  stabilizes. 

The  Interpreter’s  symbolic  methods  are  likewise  very  restricted.  It  cur¬ 
rently  uses  symbolic  algebra  only  to  derive  formulas  for  Jacobiains  and  sen¬ 
sitivities,  Another  obvious  application  for  symbolic  processing  is  to  handle 
symmetries.  For  instance,  one  should  expect  to  find  pitchfork  bifurcations 
only  if  the  system  is  invairiant  under  a  symmetry  (as  mentioned  in  section  2.1). 
M2uay  such  symmetries  can  be  identified  by  examining  the  defining  equations, 
and,  having  found  a  symmetry,  the  Interpreter  could  use  this  to  avoid  redun¬ 
dant  computations.  More  significantly,  in  the  caise  where  /  is  computed  by 
algebraic  formulas  rather  than  by  numerical  integration,  one  can  attempt  to 
recognize  and  classify  bifurcations  purely  symbolicadly,  by  deriving  the  nor- 
mad  forms  of  the  maps.  Such  computations,  using  the  Macsyma  symbolic 
adgebra  system,  have  been  demonstrated  by  Rand  and  Keith  [17]. 

It  is  in  the  au'ea  of  incorporating  more  explicit  knowledge  about  dynamics 

^®Parker  and  Chua’s  INSITE  system  [13,  14]  includes  a  wide  assortment  of  numerical 
methods  that  are  useful  in  the  study  of  nonlineu  dynamics,  including  routines  that  locate 
individual  periodic  points  and  track  them  to  bifurcations,  in  much  the  sajne  way  as  the 
tracking  phase  of  the  Interpreter  (section  3  3). 


that  the  most  important  work  remains  to  be  done.  For  instance,  beyond  iden¬ 
tifying  collections  of  fixed  points,  the  present  Interpreter  pays  no  attention 
to  the  geometry  of  the  phase  space.  Thus,  it  can  hope  to  recognize  only  those 
bifurcations  that  occur  at  individual  fixed  points.  Dealing  with  more  global 
reorganizations  of  phase  space,  such  as  saddle  connections  (see  [6])  will  re¬ 
quire  a  more  explicit  representation  of  phase-spaoe  geometry.  Additionally, 
representing  geometric  relations  between  nearby  trajectories  can  enable  a 
program  to  apply  powerful  consistency  constraints  in  guiding  its  exploration 
of  the  phase  space.  Programs  that  automatically  investigate  dynamical  sys¬ 
tems  by  “looking  at”  phase-space  geometry  have  been  demonstrated,  for 
Hamiltonian  systems  with  3  degrees  of  freedom,  by  Yip  [21,  22,  23],  and  for 
planar  vector  fields,  by  Sacks  [18]. 

A  final  important  area  for  further  work  is  to  apply  the  Interpreter  to 
systems  that  are  derived  from  models  of  physical  situations,  rather  than  pre¬ 
sented  a  priori  as  equations,  and  to  have  the  program  formulate  its  interpre¬ 
tations  in  terms  of  underlying  physical  phenomena  rather  than  only  in  terms 
of  the  bare  mathematics.  For  exaunple,  if  a  nonlinear  oscillator  describes  the 
rolling  motion  of  a  ship,  then  a  transition  to  chaos  can  be  interpreted  as  the 
possibility  of  capsizing  (see  [12]).  A  program  that  describes  chemical  oscil¬ 
lations  by  making  qualitative  investigations  of  the  corresponding  dynamical 
system  is  being  developed  by  Eisenberg  [3]. 

The  Bifurcation  Interpreter,  even  with  its  present  limitations,  lets  us 
imagine  a  new  category  of  programs  that  formulate  numerical  results  in 
qualitative  terms,  thereby  enabling  their  users  to  control  computational  ex¬ 
periments  in  terms  of  high-level  behavioral  descriptions  (see  [2]  for  further 
examples).  Over  the  past  forty  years,  increasingly  sophisticated  numerical 
tools  have  radically  transformed  the  role  of  mathematical  modelling  in  science 
and  engineering.  Looking  forward  to  combining  these  numerical  techniques 
with  symbolic  and  knowledge-based  methods,  shows  that  this  transformation 
has  only  just  begim. 
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